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Abstract 

A two-step hybrid perturbation-Galerkin method for the solution of a variety of 
differential equations type problems is found to give better results when multiple per- 
turbation expansions are employed. The method assumes that there is a parameter in 
the problem formulation and that a perturbation method can be used to construct one 
or more expansions in this parameter. An approximate solution is constructed in the 
form of a sum of perturbation coefficient functions multiplied by computed amplitudes. 
In step one, regular and/or singular perturbation methods are used to determine the 
perturbation coefficient functions. The results of step one are in the form of one or 
more expansions each expressed as a sum of perturbation coefficient functions multi- 
plied by a priori known gauge functions. In step two the classical Bubnov-Galerkin 
method uses the perturbation coefficient functions computed in step one to determine 
a set of amplitudes which replace and improve upon the gauge functions. The hybrid 
method has the potential of overcoming some of the drawbacks of the perturbation and 
Galerkin methods as applied separately, while combining some of their better features. 
In this study the proposed method is applied, with two perturbation expansions in each 
case, to a variety of model ordinary differential equations problems including: a family 
of linear two-point boundary-value problems, a nonlinear two-point boundary-value 
problem, a quantum mechanical eigenvalue problem and a nonlinear free oscillation 
problem. The results obtained from the hybrid method are compared with approxi- 
mate solutions obtained by other methods, and the applicability of the hybrid method 
to broader problem areas is discussed. 


1 Research for the first author was supported under the National Aeronautics and Space Administration 
under NASA contract No. NAS1-18107 while he was in residence at the Institute for Computer Applications 
in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665 


1 



1. INTRODUCTION 


A two-step hybrid analysis technique, which combines perturbation techniques with 
the Galerkin method, was presented and discussed by Geer and Andersen [7], where it was 
applied to some singular perturbation problems in slender body theory. That technique 
promises to be useful in the analysis of a very wide variety of differential equations type 
problems. In this paper we shall further develop the method through the use of multiple 
perturbation expansions and demonstrate its improved usefulness for several classes of 
two-point boundary-value problems for ordinary differential equations. 

The Galerkin method [6] has, of course, been known and used for a long time. But a 
principal problem associated with its successful application lies in the choice of appropriate 
basis functions. This problem was addressed when the hybrid technique was apparently 
first introduced by Ahmed K. Noor and collaborators in conjunction with the finite ele- 
ment analysis of geometrically nonlinear problems in structural mechanics. In a series of 
papers [11-21 among others] they have shown for a variety of structural mechanics and 
thermal problems that the first few terms in a Taylor series expansion of the solution of 
a parameterized system of discretized equations can be particularly effective as Galerkin 
trial functions (or basis vectors). Subsequently, Geer and Andersen [7] have shown for 
two classes of problems in slender body theory that the terms in a singular perturbation 
expansion of the solution are also effective trial functions. It has been demonstrated that 
the “reduced-basis” solutions can be useful for significantly larger values of the expansion 
parameter than the Taylor series or singular perturbation solutions on which they are 
based. A treatment of the reduced basis method from a mathematical point of view is 
given by Fink and Rheinboldt [5], 

Some general observations about the technique are the following. First, in many pertur- 
bation problems, much effort has to be expended to compute (analytically) each additional 
term in a perturbation expansion. Through the use of the proposed hybrid method, the 
known perturbation terms can be exploited more fully. Secondly, another way of viewing 
the technique is to recognize that in many perturbation expansions the functional form of 
the higher order terms can be well approximated by a linear combination of a few lower 
order terms. Thus, much of the effect of the higher-order terms may be included by using 
only a small number of lower-order perturbation functions, but with their amplitudes re- 
placed by new functions of the expansion parameter. Thirdly, since it is often possible to 
develop formal perturbation expansions for the solution about two or more values of the 
parameter, e. g. for small or large values of the parameter, the proposed method appears to 
be a convenient way of combining the information contained in these different expansions. 
The set of basis functions is chosen so as to include terms from the different perturba- 
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tion expansions. In fact, when information from two or more expansions is employed, the 
method appears to provide meaningful (and often very accurate) approximations not only 
in the regions near the expansion points but also in the “intermediate” regions for the 
parameter, e. g. in regions where the parameter is neither “small” nor “large” . 

It is our belief that the conjunction of perturbation and Galerkin techniques can be 
useful in a wide variety of application areas and in these applications the hybrid technique 
will give better approximations than the perturbation method alone. In this paper we will 
present applications of the technique to several classes of problems with second order 
ordinary differential equations. Two perturbation expansions are used in each of the 
problems illustrated. In the next section the method is described in more detail. In 
section 3 it is applied to a general class of linear two-point boundary-value problems and 
in section 4 to a related linear problem. In sections 5, 6, and 7 the method is illustrated 
by application to a representative nonlinear problem, to a quantum mechanical eigenvalue 
problem, and to a nonlinear oscillation problem, respectively. In section 8 we make some 
additional observations about the method, discuss some of the insights it provides, and 
indicate some areas for further study. 

2. DESCRIPTION OF THE METHOD 

The method we wish to describe is a two-step hybrid analysis technique. It is based 
upon the successive use of perturbation expansions and the classical Bubnov-Galerkin 
approximation technique. In the perturbation methods, an approximation to the solu- 
tion of a particular problem involving a parameter is developed in terms of a series of 
unknown functions with preassigned coefficients, i. e. gauge functions. The unknown func- 
tions are usually determined by solving a recursive set of differential equations which are, 
in general, simpler than the original governing differential equation. By contrast, in the 
Bubnov-Galerkin technique one seeks an approximate solution to the problem in the form 
of a linear combination of specified (known) coordinate (also called perturbation coordi- 
nate) functions with unknown coefficients. The coefficients are determined by demanding 
that the residual formed by substituting the trial solution into the governing differential 
equation is orthogonal to each of the coordinate functions. 

While both of these methods are useful and have been successful in providing approx- 
imate solutions to a wide variety of nonlinear (or otherwise difficult) problems, each has 
certain drawbacks. The perturbation method has at least two major drawbacks. First, as 
the number of terms in the perturbation expansion increases, the mathematical complexity 
of the equations which determine the unknown functions increases rapidly. Thus, in most 
practical applications the computed perturbation series is limited to only a few terms. A 
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second drawback to the perturbation method is the requirement of restricting the pertur- 
bation parameter to values close to the value about which the expansion was constructed 
(e. g. to small or large values of the parameter) in order to obtain solutions of acceptable 
accuracy. (These drawbacks of the perturbation method have been recognized and several 
modifications or extensions have been proposed, see e. g. Van Dyke [22] and Andersen and 
Geer [1]). The main shortcoming of the Bubnov-Galerkin method is the difficulty, from 
a practical point of view, of selecting a small number of good coordinate functions. The 
Bubnov-Galerkin method usually is not appropriate when a large number of coordinate 
functions are used because of the lack of sparseness in the generated equations. 

To illustrate the general ideas of the hybrid (or “reduced basis”) method, suppose we 
are seeking (an approximation to) the solution u to the problem 

(2.1) £(u,e)= 0, 

where £ is some differential operator and e is a parameter. (Although we shall restrict our 
attention to problems involving second-order ordinary differential equations in this paper, 
we believe that the method has a much wider range of applicability. Hence, we formulate 
the method in terms of a general differential operator £ and in section 8 discuss some of 
the other possible areas of application of the method.) Here (2.1) holds in some domain 
D , and, in addition, u must satisfy certain conditions on the boundary D. Without loss 
of generality, we can assume that these boundary conditions are homogeneous in u. 

Now the application of the hybrid perturbation-Galerkin method can be divided into the 
following two distinct steps. First, we generate the coordinate functions in a perturbation 
expansion of u about one or more specific values of the parameter e, say about e = e p , 
p = 1, 2, . . . , P. We then construct a new approximate solution consisting of a sum of some 
of these perturbation coordinate functions, each multiplied by an unknown amplitude, and 
determine these amplitudes by using the Bubnov-Galerkin method. 

To describe this idea in more detail, suppose that the solution to (2.1) can be expanded 
about c = e p into a series of the form 

(2.2) u = u * a i( c ) + °« P ( e ))> 

i=o 

where {oy(c)} is an appropriate asymptotic sequence of gauge functions and each u p - can 
be determined completely by a standard perturbation method (e. g. a composite expansion 
of inner and outer expansions) . 

A subset of all of the perturbation functions u v - are now chosen as coordinate functions 
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for the Bubnov-Galerkin technique and an approximation u for u is sought in the form 

N 

(2.3) = 

j = i 

where the (unknown) parameters 6j = Sj(e) represent the amplitudes of the coordinate 
functions u, . Here each Uj is one of the perturbation coordinate functions u p k . To determine 
the unknown amplitudes 6 ,- , we apply the Bubnov-Galerkin technique to the governing 
equation (2.1). Thus, we substitute (2.3) into (2.1) and demand that the residual be 
orthogonal to the N coordinate functions Uj over the domain D, i.e. 

(2.4) / £(%2u j 6 i ,c)u k dx = 0, k = 1,2,. . . ,N. 

Jd j=l 

Equations (2.4) represent a set of equations for the unknown amplitudes 6j. While (2.4) 
must, in general, be solved numerically, doing so is much simpler than numerically solving 
(2.1). In particular, for a fixed value of e, the solution to (2.4) is a point in iV-dimensional 
space, where N is reasonably small, while the solution of (2.1) is a continuous function. 

We should note that this particular choice of coordinate functions overcomes the main 
drawback of the Bubnov-Galerkin method. By the way they are constructed, the pertur- 
bation coordinate functions are (under certain assumptions) elements of a set of functions 
which span the space of solutions in a neighborhood of their point of generation. Thus, 
they should fully characterize the solution u in that neighborhood. Also, in many appli- 
cations, the functions u p are determined by solving a set of linear equations, even though 
the original operator £ may be nonlinear. The first property (completeness) is necessary 
for the convergence of the Bubnov-Galerkin method, while the second property (linearity) 
enhances the effectiveness of the proposed hybrid method for solving nonlinear problems. 

Another important property of the proposed method is that the coordinate functions, 
i.e. the perturbation functions, do not need to come from a regular perturbation expansion. 
In fact, all that is needed is a formal asymptotic expansion of the solution to (2.1) for 
values of £ near e p in the form of (2.2), where {o^(e)} is a sequence of appropriate gauge 
functions, e. g. expressions which may involve log(e) and/or fractional powers of e. Thus, 
the proposed method may be applied to singular as well as regular perturbation problems. 

Still another important property of the proposed method is its ability to combine 
information about the solution from several different asymptotic expansions. As we shall 
see in the following sections, the hybrid solution, at least in all of the examples we examined, 
has the nice property of providing reasonable (and often, very accurate) approximations 
to the solution for all values of the parameter which lie between the points about which 
the formal asymptotic expansions were constructed. 
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We shall now illustrate the hybrid method by applying it to several different types of 
model two-point boundary-value problems with second-order equations. 


3. LINEAR EQUATIONS 

In this section, we shall apply the hybrid method to a general linear second order 
problem of the form 

(3.1) e u" + a(x) u' + b(x) u = c(x) , 0 < x < 1, 

u(0) = 0 = u(l). 

For simplicity, we shall assume that a(x) > 0 on 0 < x < 1, so that, for small positive 
values of e, the boundary layer will form near x = 0. 

Using standard perturbation techniques (see, e. g. [10]), it is possible to construct 
formal asymptotic expansions of the solution to problem (3.1) for both small or large values 
of e. Thus, following the notation of the previous section, we shall consider expansions of 
u about P = 2 values of e, namely, £i = 0 and £2 = 00 . In particular, for small values of e, 
the expansion is of the form (2.2) with aj = c J , and 


(3.2) u) = yj {x) + Yj{x/e) - £ y/l^O) (*/ c)*'/*! - * 


1=0 




i=0 


The first term on the right side of (3.2) is the j th term in the outer expansion of (3.1), 
while the second term is the j th term in the inner expansion. The next term is the j tk term 
in the outer expansion of the inner solution. Thus, the first three terms together represent 
the j th term in a uniform (composite) expansion of the solution for small values of e. This 
j th term vanishes at x = 0, but does not vanish at x = 1, where it is exponentially small. 
Thus, the final term, which is exponentially small when e is small, is added to insure that 
u l j vanishes identically at both 0 and 1. In particular, we find 


(3.3) 


»(*) = ffwwy* 

E(x) = 


y„(i/e) = y 0 (o){i-e-*( 0 )*/‘}. 

Other terms can be computed in a straightforward manner (see, e. g., [10]). 

For large values of e, the expansion of the solution of (3.1) is again of the form (2.2) 
with a 2 j = £ _,_1 , and 


(3.4) 


u 2 — Zj(x) — Fj(x) xFj( l). 
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y>i. 


F A X ) = f \ fj{t)dtds , 

JO Jo 

f 0 {x) = c(x), fj(x) = —a(x) Zj_i(x) - b(x) zy_i(x), 

We now use these two expansions to construct a hybrid solution to (3.1-). In particular, we 
select n i terms from the small-e expansion (3.2) and n 2 terms from the large-e expansion 

(3.4) and write our hybrid solution in the form (2.3), with N = ni + n 2 . We then substitute 
this expression into (3.1) and use the orthogonality conditions (2.4) to obtain the conditions 

N 

(3.5) £4y(0*i( e ) = B k, k = l,2,...,N 

;= i 

At,- = J (e u'j + a u'j + b uy) u* dx, Bk = J cut dx. 

Thus, we see that the 6^ s satisfy a set of linear algebraic equations, which, of course, follows 
from the fact the original problem (3.1) is linear. Once these equations have been solved, 
the values of the 6j' s can be used along with the corresponding (coordinate) functions Uj 
in (2.3) to yield our hybrid solution. 

To illustrate our method for this class of problems, we apply it to the following three 
model problems. The first problem is 

(3.6) c u" + u' + u = sin(27r x), 

u(0) = 0 = u(l). 

In Figure 1, we have plotted for e = 0.2 in equation (3.6), (a) the singular perturbation 
solution P 0 [l] based on the one-term expansion about e = 0, (b) the regular perturbation 
solution Poo[l] based on the one-term expansion about e = oo, (c) the hybrid solution H[l,l] 
which unifies the two perturbation expansions, and (d) the “exact” solution computed by 
a numerical shooting method. (Here the notation fT[ni,n 2 ] stands for a hybrid solution 
based on nj terms from the first perturbation expansion and n 2 terms from the second 
expansion.) While neither one-term perturbation expansion succeeds well near e = 0.2, 
the hybrid solution based on just two perturbation coefficient functions gives a reasonable 
approximation not only at e — 0.2 but, in fact, for all positive values of e. 

The second problem is 

(3.7) £u" + «" I u' + e I u = l, 

u(0) = 0 = u(l). 

In Figure 2, we have plotted for e = 0.2 in equation (3.7) comparable curves to those 
shown in Figure 1. Again while neither one-term perturbation expansion succeeds well 
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near e = 0.2 the hybrid solution based on just two perturbation coefficient functions again 
gives a reasonable approximation for all positive values of e. 

To illustrate a problem wherein the H[ 1, 1] solution is not adequate we consider the 
equation 

(3.8) eu" + e x u' + u — 1 + t~ z sin(27rx), 

u(0) = 0 = u(l). 

In Figure 3, we have plotted for c = 0.3 in equation (3.8) the following: (a) the singular 
perturbation solutions Po[l] and Po[2] based on the one-term and two-term expansions, 
respectively, about c = 0, (b) the regular perturbation solutions P<»[1] and P 00 [2] based 
on one-term and two-term expansions about e = oo, (c) the hybrid solutions H[ 1, 1] and 
H[ 2,2], and (d) the “exact” solution. Here H[ 1, 1] gives a poor approximation for e = 0.3. 
However the use of one additional perturbation coefficient function from each expansion 
results in a reasonable approximation at e = 0.3 and in fact for all positive values of e. We 
shall return to these examples in section 8. 

4. A SPECIAL SINGULAR LINEAR PROBLEM 

The hybrid method can be applied successfully to a wider range of linear problems 
than considered in the previous section. To illustrate this point we shall apply the method 
to the model problem 

(4.1) eu" + xu' — xu = x(a + bx), 0 < x < 1, 

u(0) = 0 = u(l), 

where a and b are specified constants. (See Hanks [8] where this equation is used to model 
heat flow near an ocean rise.) In (4.1) the coefficient of u' vanishes at x = 0. Thus, even 
though a boundary layer still forms near x = 0, the formal perturbation solution for small 
positive values of e contains fractional powers of e as well as log(e). 

For small values of e it is straightforward to show that 

uj(x, e) + 0(e 1/2 ) -I- 0(e log(e)) as e -> 0, 
u 0 (x,c) - xu 0 (l,c), 

-bx + (^)(e* - 1) + cr^dt 

^{ 212 * - 0 - 6 } . 


u{x , 6) = 
uj(x,€) - 

(4.2) 

&o(x,€) = 

c = 
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In (4.2) u 0 is the first term in a composite expansion based upon a one-term inner ex- 
pansion and a one-term outer expansion. Then uj is obtained by subtracting from u 0 the 
(exponentially small) term xuo(l,e) to ensure that the boundary conditions are satisfied 
exactly. 

For large values of e, u has a perturbation expansion in the form of (2.2) with a* = 

The coefficients uj are again given by (3.4) with a(x ) = x, b(x) = -x, and c(z) = x ( a+bx ). 

We can now apply our hybrid method in the manner outlined in the previous section 
and obtain (linear) equations of the form (3.5) to determine the coefficients Sj. In Figure 

4, we have plotted (for e = 0.3) the approximate solutions to (4.1) based upon a one- 
term small-e perturbation solution (from equations (4.2)), a one-term large-e perturbation 
solution (from equations (3.4)), the hybrid solution H[l,l] based on these two one-term 
expansions, and a solution obtained by a numerical shooting method. Clearly, the hybrid 
solution is more accurate than either perturbation solution alone. In fact, the value of 
e used in this figure is very near the value of e corresponding to the poorest agreement 
between the hybrid and numerical solutions (in the relative L 2 sense) for all positive values 
of e. We shall discuss this example further in section 8. 

5. A NONLINEAR EXAMPLE 

To illustrate the application of the hybrid method to a nonlinear problem, we consider 
the problem of determining the excess pressure u(x) under a slider bearing, which was 
discussed for large bearing numbers by DiPrima [4]. u(x) satisfies 

(5-1) c^|fi 3 (x)[l + u(x)]^| = ^{/i(x)[l + u(x)]}, 0 < x < 1, 

u(0) = 0 = u(l). 

In (5.1) e is the inverse of the bearing number and h(x) is the (prescribed) film thickness 
under the bearing, which has been normalized so that /i(l) = 1. For small values of e 
DiPrima constructed an asymptotic expansion of the form (2.2) with e x — 0 and aj = e>, 
j = 0, 1, In particular, using his results, we find 


(5.2) 


+ Q - MO) - (i - *) [<2 Q) - MO) + 1] , 


h{0) log 


h(0) - 1 -Q 
h{0) - 1 


+ Q — 



The first three terms in the expression for uj represent the first term in a (uniform) 
composite expansion of the solution, while the last term, which is exponentially small for 
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small values of e is added to ensure that uj satisfies the boundary conditions exactly. The 
second expression in (5.2) defines Q implicitly as a function of (1 — x)/e. 

For large values of e, u has a regular perturbation expansion of the form (2.2) with 
a 2 = In particular, we find 

(5.3) Uq = [ h~ 2 (t)dt + C f h~ s (t)dt , 

Jo Jo 

C = — J h~ 2 (t)dt j J h~ s (t)dt. 

More terms in this expansion can be computed in a straightforward manner using simple 
quadratures. 

We now select n\ terms from the small-e expansion of u and n 2 terms from the large-e 
expansion and form our hybrid solution u in the form (2.3), where N = n\ + n 2 . The 
equations (2.4) for the amplitudes {£,-} are now quadratically nonlinear. However, we note 
that, for small values of e, the fy’s corresponding to small-e perturbation coordinate func- 
tions should be approximately e J , while the 6/s corresponding to the large-e perturbation 
coordinate functions should be approximately zero. Thus, these equations can be solved 
easily (and efficiently) by using Newton’s method, starting at small values of e and then 
proceeding to larger values of e. 

To illustrate the accuracy of our hybrid method, we apply it first to the case of a 
linear film thickness (i. e., a wedge slider) described by h(x) = 1 + a(l — x). In Figure 5, 
we have plotted the pressure distribution obtained from (5.1) by using a one-term small- 
e perturbation expansion P 0 [l], a one-term large-e expansion P^l], the hybrid solution 
H[ 1, 1] based on the two coordinate functions contained in these one-term expansions, 
and a solution obtained by a numerical shooting method. As can be seen, the agreement 
between the hybrid solution and the numerical solution is excellent, and the hybrid solution 
is considerably better than either the small-e or large-e perturbation solutions alone. In 
fact, as we performed similar comparisons with this example for other values of e, we 
found that the value of e = 0.7, corresponding to Figure 5, represented the “poorest” 
agreement between the hybrid and numerical solutions, with a relative L 2 -e rror of about 
1%. In Figure 6, we have plotted the hybrid solution jCT[ 1, l] and the numerical solution 
for several different values of e. For values of c substantially larger or smaller than those 
shown in the figure, the relative L 2 -error is considerably less than 1 %. In Figure 7, we 
have plotted <5i and <5 2 (the coefficients of the small- and large-e coordinate functions, 
respectively) as functions of e. As the figure illustrates, Si approaches 1 and 82 approaches 
0 as e approaches zero, while S 2 dominates S\ as e approaches infinity. 

In Figure 8, we have made a similar comparison to that shown in Figure 5, except 
that the film thickness corresponds to a parabolic shape as discussed by DiPrima. The 
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results depicted in Figure 8 are typical of those we obtained for several different parabolic 
shapes. For the shape chosen, the value of c used in the figure corresponds to the “poorest” 
agreement between the hybrid and numerical solutions. 

In section 8, we shall discuss some of the insights this example provides for the appli- 
cation of the hybrid method to other nonlinear problems. 

6. AN EIGENVALUE PROBLEM 

We feel that our hybrid method can be particularly useful when applied to eigenvalue 
problems and we are currently investigating its application to some general classes of such 
problems. However, for our present purposes, we shall be content with illustrating its 
application to second order eigenvalue problems by considering the following (quantum 
mechanical type) model problem: 

£(u, e, A) = u"(x) + {A - V (x, c)} u(x) = 0, -oo < x < oo, 

/ oo 

u 2 (x) dx — 1, 

-oo 

V (x, e) = x 2 + e x 4 . 

Here A is an eigenvalue and V is a potential function. When e = 0, (6.1) reduces to a 
well-known quantum mechanical harmonic oscillator problem with a quadratic potential 
V(x, 0) = x 2 . For positive values of e we expect to find eigenvalues and corresponding 
eigenfunctions which reduce to those of the harmonic oscillator as e goes to zero. However, 
for negative values of c, there are no solutions to the problem (6.1), because V becomes 
negative for large values of x. Consequently, we expect that the classical perturbation 
solution to (6.1) will have a zero radius of convergence. This is indeed the case, as we shall 
now illustrate. 

For small values of e, the classical perturbation solution to (6.1) is of the form (2.2), 

i.e., 

(6.2) u=Yl u ; (x) e y + 0(e N ), A = Y, X i e ’ + °( e ")• 

i = o i = o 

Substituting (6.2) into (6.1), we find 

(6.3) A 0 = 2fc + 1, u 0 (x) = e~ x3 / 2 hk{x), k = 0, 1,2, . . . , 

= / fi( x ) «o(z) dx, Uj{x) = pj(x) e" z2/2 , 

J — OO 
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i-1 

fj(x) = x 4 tlj—i(x) - “/-<(*)» / > 1- 

»•= 1 

In (6.3), h*(x) is the Hermite polynomial of order k and p, is a polynomial of order 
Aj + k. For our purposes here, we shall only consider perturbations about the ground state 
eigenvalue A = 1 corresponding to k — 0. Similar results hold for the higher eigenvalues and 
some of these results are reported and discussed in [2]. In particular, using the symbolic 
manipulation system MACSYMA [9], we find (for A; = 0) 

(6.4) A = 1 + (3/4)e - (21/16) e 2 + (333/64)e 8 - (30885/1024) e 4 , 

h 0 (x) = 1 /tt 1 / 4 , Pi(x) = — (4x 4 + 12x 2 — 9)/(327r 1 / 4 ), 

p 2 (x) = (48x 8 + 416x 6 + 1272x 4 + 3384x 2 - 4677)/(6144?r 1 / 4 ). 

(Several more terms in this series have been computed using a simple BASIC program. 
Both the ratio and root tests clearly indicate a zero radius of convergence for the pertur- 
bation expansion.) 

We now apply our hybrid method to this problem using at first only terms from the 
small-e perturbation solution (6.2). Thus we look for an approximate solution for u of the 
form (2.3), where, for convenience, we now begin the summation at j = 0 and let the upper 
limit be N — 1. Then each Uj is defined in (6.3) with k = 0. Substituting this expression 
into (6.1), the orthogonality conditions (2.4) can be written as the matrix-vector eigenvalue 
problem 

(6.5) (C + A£>)£=0, 

where 6 is the vector of the unknown amplitudes 6j and the N by N matrices C and D 
are defined by 

Cij = J {«;'(x)-(x 2 + 6x 4 )u J (x)}u,(x)dx, 

(6.6) 

/ oo 

u,(x) Uj(x) dx, i,j = 0, 1, . . . , N — 1. 

-OO 

Thus, for a fixed value of e, the smallest eigenvalue of (6.5) will provide our approximation 
to the smallest eigenvalue of (6.1), while the corresponding eigenvector, when used in (2.3), 
will provide our approximation to the corresponding eigenfunction. 

In Figure 9, we have plotted the approximations to the ground state eigenvalue for 
different values of c obtained by using N = 1, 2, . . . , 6 terms in the perturbation and hybrid 
solutions, along with the eigenvalues obtained by a purely numerical solution to (6.1). The 
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figure indicates that the perturbation solution is essentially useless for approximating the 
eigenvalues, except for extremely small values of e. It also indicates that for any fixed 
value of e the hybrid approximations appear to be converging to the correct eigenvalue 
as N increases, although, as expected, the convergence is slower for larger values of e. 
Incidentally, the two-term perturbation solution is identical to the hybrid solution with 
N = 1. 

To find approximations to the solution of (6.1) valid for large values of e, it is convenient 
to introduce the new variables n = e -1 / 3 A, z = e 1 / 6 x, and v(z ) = e 1 / 12 u(x). Then 
(6.1) becomes 

(6.7) v"(z ) + j/x — z* - e _2 / 3 z 2 1 v(z ) = 0, 



Then v(z, e) has a perturbation solution in the form of (2.2) with a 2 (e) = e -2j / 3 , with a 
similar expansion holding for /x = /x(e). In particular, the leading term in this expansion 
satisfies the problem 

v'q[z) + {/x 0 - z 4 } v 0 (z) = 0 

/ OO 

vl(z)dz = 1. 

-OO 

Although we have not been able to express the solution to (6.8) in closed form, it can be 
solved numerically in a straightforward manner. In particular, we find that /x 0 = 1.0604, 
which implies that A = 1.0604 c 1 / 3 + 0(e -1 / 3 ) as e — ► oo. 

We now add the function vo{z) to the set of basis functions formed in (6.3) and construct 
some new hybrid approximations. That is, we use v 0 (z) as well as some of the small- 
e perturbation coordinate functions in our approximation (2.3). In Figure 10 we have 
plotted the eigenvalues obtained from the one-term large-e perturbation solution P<»[1] 
as well as those obtained from the hybrid solutions H[0, l], H [1, 0 ], and #[1,1]. As the 
figure clearly indicates, the inclusion of one term from the large-e perturbation expansion 
dramatically improves the accuracy of our approximations. 

The quality of the hybrid approximations to the corresponding ground state eigenfunc- 
tion is illustrated in Figure 11. Here the hybrid approximations based on different numbers 
of terms from the small-e and large-e perturbation solutions are plotted along with the 
eigenfunction obtained by a numerical shooting method for a fixed value of e. Again, the 
hybrid approximations based on just the small-e perturbation solution, i.e., H[M, 0], are 
apparently converging as M increases, while the inclusion of one term from the large-e 
expansion again dramatically increases the accuracy of our hybrid approximation. 
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7. A NATURAL FREQUENCY CALCULATION 

The hybrid method we are proposing can also be useful in determining approximations 
to the period T (or, equivalently, the natural frequency 2n/T) for a nonlinear oscillator 
involving a parameter e. To illustrate applications to problems in this area, we consider 
\ the problem of determining the natural frequency of a simple mechanical system consisting 

of a mass m restrained by two identical springs, each having natural length L and spring 
^ constant k, midway between two parallel walls a distance 2 d apart, as discussed by Arnold 

and Case [3]. Thus, if we let h be the maximum displacement of the mass, the displacement 
I along the centerline between the two planes is described by hu(x), where 

(7.1) i/ 2 u" + u + fxu |l — = 0, 

u = 1, u' = 0 at x = 0; u(x + 2n) = u(x). 

In (7.1), e — h/d, n = A/(l — A), A = L/d < 1, v — 2n/Tw, w = yj2k(\ — \)/m and 
x = wut, where t is time. 

For small values of e, u has an expansion of the form (2.2), with a) = e 2j , and a 
similar expansion holds for u 2 . In particular, using the Lindstedt-Poincare method (see, 
e. g. Nayfeh [10]), we find 

i 

| (7.2) u = cos(x) + (ji/64) {cos(3x) — cos(x)} t 2 + 0(e 4 ) 

v 2 = 1 + (3/»/8)e 2 + 0(e 4 ) as £-*0. 

I For large values of e, u has a (composite) expansion in the general form of (2.2) with 

a 2 = e~\ in which the formal outer expansion, valid for cu > 1, must be supplemented by 
inner expansions around x = tt/2 and Ztt/2, where u vanishes. In this case, it is easy to 
show that u is still well approximated by the first term on the right side of (7.2), while 

(7.3) i/ 2 = 1 + n + 0(e~ 1 ) as e — ► oo. 

• To apply our hybrid method to this problem, we look for an approximate solution u 

for u in the form 

l N 

(7.4) u = Uj(x) Sj, S 0 = 1, 

! j=o 

| where uo(x) = cos(x) and the remaining uy(x) can be selected from either the small-e or 

large-e expansions of u. Substituting (7.4) into (7.1), we see that the initial conditions and 
the periodicity condition are satisfied, while the orthogonality conditions (2.4) become 

1 (7.5) F j (v 2 ,6 1 ,6 2 ,...,6 n ) = 0, j =0,1,..., AT, 
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u 2 u" + u + /xu( 1 — 


\/TT 


=p)}«J ds - 


Equations (7.5) are a system of N + 1 equations for the N + 1 unknowns i/ 2 and 6j, 
j = 1 , 2 ,..., AT, which can be solved using Newton’s method. For the .special case when 
N = 0, we can solve (7.5) explicity and find 


(7.6) 



+ € 2 COS 2 (x) 


cos 2 (x) dx. 


If we can expand (7.6) for small values of e, we recover the expansion for u 2 in (7.2), while 
if we expand it for large values of e we recover the expansion (7.3). In Figure 12, we have 
plotted v 2 determined by (7.6) as a function of e, and have also plotted selected values 
determined by purely numerical means. As the figure indicates, the agreement between 
the hybrid and numerical results is excellent. 


8. OBSERVATIONS AND CONCLUDING REMARKS 

Sections 3 through 7 illustrate that the hybrid perturbation-Galerkin method unifying 
two perturbation expansions can successfully be applied to a wide variety of two-point 
boundary-value problems. 

In Section 3 a whole class of linear problems is treated. Three examples are chosen from 
this class for purposes of illustration. In each problem a singular perturbation expansion 
is constructed about 6 = 0, a regular expansion is constructed for large values of e, and 
the hybrid method combines terms from both expansions. For many problems in this 
class, as exemplified by equation (3.6) and Figure 1, or equation (3.7) and Figure 2, 
the #[ 1, 1] type solutions provide reasonable approximations for all positive values of e. 
Probl ems where the #[l, 1] approximation performs poorly seem to be somewhat unusual. 
Equation (3.6) and Figure 3 represent just such a case. For e = 0.3, Figure 13 compares 
hybrid solutions given by #[2,0], #[1,1], and #[0,2] and an “exact” solution determined 
numerically. In this case #[0,2] is the best of the three hybrid solutions. This is fortunate, 
since the regular perturbation functions used in #[0, 2] are much easier to compute than 
the singular perturbation functions needed in #[2,0]. However, it is unfortunate in that 
we did not yet know how to determine a priori which of these approximations will have 
better accuracy. 

To compare the accuracy of various perturbation and hybrid solutions the correspond- 
ing relative L 2 -error norms are shown in Figure 14 as functions of e. The typical behavior 
is 
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(i) the hybrid solution #[#,0] for each N = 0,1,2,... is more accurate than the 
perturbation solution Po[N] upon which it is based; 

(ii) both #[#,0] and Po[N] tend to be more accurate for small e than for large e; 

(iii) the hybrid solution #[0, iV] for each N is more accurate than the perturbation 
solution #<»[#] upon which it is based; 

(iv) both #[0,iV] and Poo[^] tend to be more accurate for large e than for small e; and 

(v) hybrid solutions based on two expansions, such as H[N,N\, tend to be accurate 
everywhere for sufficiently large N. 

The problem (3.6) seems to be somewhat atypical. For instance, it is observed that 
Poo[l] outperforms #[0, 1] and even #[1, l] for some values of e. However, #[ 0,2] is much 
better than Poo[2] and #[2,2] is much better that #[0, 2] or #[2,0]. 

In Section 4 we examined a linear problem whose small-e perturbation expansion in- 
volves fractional powers of e as well as terms involving log (e). For this example, the results 
of the hybrid method were still very good, with the hybrid solution based upon just one 
term from each of the two perturbation expansions providing a good approximation to the 
exact solution. This example provides a deeper insight into our method in the following 
way. If we consider a hybrid solution based on only one term from the small-e expan- 
sion, i. e. u = 5iuJ(x,e), where uj(x, e) is defined in (4.2), then from (3.5) we see that 
= Bi/Aiu were B\ and An are defined in (3.5). Using the definition of uj(x, e) in this 
expression, we find that has an expansion for small e in the form of a series in powers 
of e 1 / 2 . Thus, some of the higher order e dependence of the solution is “anticipated” by $i. 
Moreover, if we consider a hybrid solution based on two terms from the small-e expansion, 
i. e. if we include the coefficient uj(x, e) of e 1 / 2 in (4.2) in our set of basis functions, we find 
that both Si and S 2 have expansions for small e that not only involve powers of e 1 / 2 but also 
involve log(e). This last fact stems from the fact that uj(x,e) involves an integral whose 
integrand behaves like x -1 near x = 0 and whose lower limit of integration is proportional 
to e 1 / 2 . Thus, again, some of the higher order singular behavior of the solution seems to 
be accounted for in the hybrid coefficients. 

Since the class of nonlinear differential equations is too broad for a general treatment, a 
specific nonlinear differential equation was selected for treatment in Section 5. The hybrid 
method again combines a singular perturbation expansion about e = 0 with a regular 
expansion valid for large e. Very good accuracy is obtained by using only one term from 
each of the two expansions, but, of course, better accuracy is expected if more terms are 
used. As is generally the case for nonlinear problems, the equations which determine the 
Sj are nonlinear and require a numerical solution. However, these equations can be solved 
efficiently using Newton’s method. For small values of e the form of the Sj is known. Thus, 
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there is a good starting point for e near zero and one can incrementally increase e, solving 
for the Sj one step at a time. 

An eigenvalue problem arising in quantum mechanics is discussed in Section 6. For this 
problem the well known perturbation expansion at e = 0 of the ground state wave function 
consists of terms which are polynomials in x multiplied by exp(— x 2 /2). However, this 
expansion has zero radius of convergence, and perturbation expansions of order greater 
than one are essentially useless. On the other hand, the hybrid solutions based on a few 
terms in this divergent expansion (see Figure 9) give accurate results for e from 0 up to 
order of magnitude one. While the hybrid solutions appear to converge for all e, the rate 
of convergence is slower as e increases, so the technique, while a decided improvement over 
the perturbation method, is, in fact, useless for high values of c. In Andersen and Geer 
[2] several low energy states were examined with this method and the same conclusions 
apply except that the convergence is observed to be slower for higher energy levels. A 
perturbation expansion for large values of e can also be determined even though none of 
its terms can be expressed explicitly in terms of elementary functions. However, they can 
be determined numerically. Hybrid solutions based on just the first term of the expansion 
at e = oo and on one or more terms of the expansion at e = 0 give results which are far 
more accurate than the hybrid solutions based solely on terms from the expansion about 
c = 0 (see Figures 9 and 10). The use of the hybrid method with the inclusion of one term 
from the large-e expansion gives useful approximations for all positive values of e. 

Finally, in Section 7 there is a natural frequency calculation for a mechanical system 
consisting of a mass and two springs. The perturbation expansions are performed for small 
and large amplitudes of the oscillations. Again by using only one term from the small- 
amplitude expansion and one term from the large-amplitude expansion, a hybrid solution 
is constructed which yields better approximations for the natural frequency than do the 
perturbation solutions on which it is based. 

Although we do not as yet have a general proof that our hybrid method will converge to 
the exact solution as the total number of terms N in the solution increases, we are currently 
working on such a proof for general classes of linear problems, such as those considered in 
section 3. As a preliminary step in this direction, a simple linear example with constant 
coefficients was analyzed in depth in Andersen and Geer [2] and a general expression was 
determined for all of the S } - associated with an arbitrary value of N. Examples such as 
these, where the details can be closely examined, may provide insight as to how general 
convergence proofs may in the future be constructed. We also leave consideration of rates 
and types of convergence to future investigations. 

The hybrid method in any previous publications known to us has been based on pertur- 
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bation expansions about a single value of e. In this paper we show several examples where 
two perturbation expansions, based on widely separated values of e, are used. Work in- 
volving expansions about three or more points is in progress and will be reported in future 
papers. On the basis of the examples studied here, it seems that when expansions about 
two points are feasible, better results are obtained for values of the parameter between 
these two points when terms from both expansions are used rather than when terms from 
only one expansion are used. 

The hybrid method may be generalized in the sense that the Bubnov-Galerkin step, 
wherein the same set of functions is used as “trial” and “test” functions, can be modified 
as to use different trial and test functions, or even other weighted residual techniques 
including collocation. Such generalizations are planned for future work. Still other future 
work will include systems of equations. 

The hybrid method based on single perturbation expansion is being applied to a very 
broad class of problems which includes semi-analytic solutions not only to ordinary differ- 
ential equations with one parameter, but also to integral equations, to partial differential 
type equations, to expansions in more than one parameter, to problems discretized by 
finite differences or finite elements, and to problems with broken symmetry. The research 
presented here suggests that the hybrid method using multiple perturbation expansions 
may be useful in these areas as well. 
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FIGURE CAPTIONS 


1. Perturbation (P 0 [l] and Poo[l]), hybrid (#[1,1]) and exact (circles) solutions for the 
two-point boundary- value problem (3.6) with e = 0.3. 

2. Perturbation (Po[l] and Poo[l]), hybrid (#[1,1]) and exact (circles) solutions for the 
two-point boundary- value problem (3.7) with e — 0.2. 

3. Perturbation (Po[l]j P > o[2], Poo[l] and P<»[2]), hybrid (#[1,1] and #[2,2]) and exact 
(circles) solutions for the two-point boundary-value problem (3.8) with e = 0.3. 

4. Perturbation (Po[l] and Poo[l]), hybrid (#[1,1]) and exact (circles) solutions for the 
two-point boundary- value problem (4.1) with a = b = 1 and e = 0.3. 

5. Perturbation (Po[l] and Poo[l])j hybrid (#[1,1]) and exact (circles) solutions for the 
wedge slider bearing problem (5.1) with a = 1 and e = 0.7. 

6. Hybrid (#[1,1]) and exact (circles) solutions for the wedge slider bearing problem 

(5.1) with a = 1 and several values of e. 

7. The amplitudes and 62 in the hybrid solution #[1,1] for the wedge slider bearing 
problem (5.1) with a = 1. 

8. Perturbation (Po[l] and #<*,[1]), hybrid (#[1,1]) and exact (circles) solutions for the 
slider bearing problem (5.1) with a parabolic film thickness, h(x) = 1+2(1 — x) 2 and 
e = 0.5. 

9. Regular perturbation (P 0 [l] through Po[6]), hybrid (#[1,0] through #[6, 0]) and 
exact (circles) solutions for the lowest eigenvalue of the quantum mechanical problem 

( 6 . 1 ) . 

10. Perturbation (Poo[l]), hybrid (#[0,1], #[2,0] and #[1,1]) and exact (circles) solu- 
tions for the lowest eigenvalue of the quantum mechanical problem (6.1). 

11. Hybrid solutions (#[1,0] through #[4,0] and #[l,l]) and exact (circles) solutions 
the eigenfunction corresponding to the lowest eigenvalue of the quantum mechanical 
problem (6.1) with e = 0.3. 

12. Square of the natural frequency in the nonlinear oscillator problem (7.1) for several 
values of the parameter /z as computed by the hybrid method (equation (7.6)) and a 
numerical method. 
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13. Two-term hybrid (#[0,2], H [1, l] and #[2,0]) and exact (circles) solutions to equa- 
tion (3.6) for 6 = 0.3. 

14. A comparison of relative L 2 -error norms as functions of c for several perturbation 
and hybrid solutions to equation (3.6). 
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